Nonlinear dynamics for vortex lattice formation in a rotating Bose-Einstein 

condensate 



Kenichi Kasamatsu and Makoto Tsubota 
Department of Physics, Osaka City University, Sumiyoshi-Ku, Osaka 558-8585, Japan 

Masahito Ueda 

Department of Physics, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8551, Japan 

(Dated: February 1, 2008) 

We study the response of a trapped Bose-Einstein condensate to a sudden turn-on of a rotating 
drive by solving the two-dimensional Gross-Pitaevskii equation. A weakly anisotropic rotating po- 
tential excites a quadrupole shape oscillation and its time evolution is analyzed by the quasiparticle 
projection method. A simple recurrence oscillation of surface mode populations is broken in the 
quadrupole resonance region that depends on the trap anisotropy, causing stochastization of the 
^vq . dynamics. In the presence of the phenomenological dissipation, an initially irrotational condensate 

is found to undergo damped elliptic deformation followed by unstable surface ripple excitations, 
some of which develop into quantized vortices that eventually form a lattice. Recent experimental 
results on the vortex nucleation should be explained not only by the dynamical instability but also 
by the Landau instability; the latter is necessary for the vortices to penetrate into the condensate. 
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^ . I. INTRODUCTION 

Since realization of a Bose-Einstein condensate (BEC) of alkali-metal atomic gases, much attention has been focused 
on the dynamical phenomena associated with superfluidity. The remarkable feature reflecting superfluidity appears 
in the response to external rotation. Recent observation of a quantized vortex lattice in trapped BECs ||, ]ij 

■ confirmed the evidence of superfluidity. Madison et al. observed directly the nonlinear dynamical phenomena such 
as the vortex nucleation and lattice formation in the rotating condensate ||. Such visualized results have greatly 

■ contributed to elucidation of static and dynamic properties of quantized vortices. 

The dynamics of dilute BECs has been successfully described by the Gross-Pitaevskii (GP) mean field model. For 
the quantized vortices in a trapped BEC, various theoretical studies have been made based on this model ||. The 
mechanism of vortex nucleation in rotating trapped BECs is one of the important topics. Vortex nucleation of this 
system differs from that of a superfluid helium system in the ratio of the coherence length £ to the system size L. In 
the former where £ < L, vortex nucleation is related to the instability of collective excitations whose energy scale is set 
\ by the confining potential. In the latter where £ <C L, it is determined by the local dynamics. A number of theoretical 
papers have discussed possible mechanisms of vortex nucleation @, |, [l2[|l3[0[l5l00^0|^|l). 



However, only a few papers made full numerical analysis of the time-dependent GP equation, which is necessary to 
understand the results in Re£[ 5| . Although the imaginary time propagation of the GP equation is a powerful scheme 
to find equilibrium states P, |l~3|, the dynamical process toward such states cannot be revealed by this method. Feder 
et al. solved numerically the time-dependent GP equation in a rotating frame, but the motion of generated vortices 
remains turbulent, forming no vortex lattice [ ^0[ . Tsubota et al. p"4[ ] included a phenomenological dissipation into 
the GP equation to simulate a vortex lattice formation, obtaining the results consistent with those of Madison et al 
The whole story of the vortex lattice formation has been clarified as follows: (i) the condensate makes elliptic 
deformation to oscillate, (ii) surface waves are excited at the boundary of the condensate, (iii) quantized vortices enter 
the condensate from the boundary, forming a lattice. 

It has been not yet clear what relation between the dynamical processes (i)-(iii) and the intrinsic instability of a 
rotating condensate. There are two considerable instabilities for vortex nucleation, namely, the dynamical instability 
and the Landau instability @, §, [l5[ @ |l|, The former originates from the imaginary frequency of the 
excitation mode, giving rise to the exponential growth of the unstable mode even in the energy-conserving dynamics. 
The latter occurs when the corresponding excitation frequency becomes negative in the reference frame and some 
dissipation works. These two instabilities work generally in different parameter ranges. Thus, one may question 
which instability is important for the actual experiments on vortex nucleation. It should be noted that the vortex 
nucleation frequency in a series of ENS experiments |l], || coincides with the prediction based on the dynamical 
instability ]l2| , contrary to some analyses based on the Landau instability. 

In this paper, we investigate theoretically the detailed dynamics of a BEC subject to external rotation by the 
numerical analysis of the two-dimensional GP equation, and solve the above questions. The first issue of this paper 
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deals with a response of a BEC to a sudden turn-on of a rotating drive within the energy-conserving dynamics. 
The rotating potential excites chiefly the quadrupole surface mode with angular momentum I — 2 and distorts the 
condensate into an ellipse. Because of the anisotropy of the trapping potential and the nonlinear atomic interaction, 
the different surface modes of oscillation are coupled to each other, causing complicated nonlinear dynamics. To clarify 
the mode coupling we use the quasiparticle projection method E2 ], which allows us to decompose the macroscopic 
wave function into the condensate and noncondensate parts and determine the populations of each mode. We find 
that the condensate makes the simple periodic oscillation for most values of the rotation frequencies; populations in 
excited modes are recovered completely to initial values in the sense of the Fermi-Pasta- Ulam recurrence |23f| . In 
resonance of the quadrupole mode, however, the simple recurrence is broken into chaotic behavior, reflecting the 
dynamical instability of the rotating condensate |L2| . The increase of the trap anisotropy extends the range of the 
rotation frequency for the resonance excitation. The chaotic dynamics yields violent density and phase fluctuations 
at the condensate surface. The generated surface ripples slightly increase of the total angular momentum, but they 
never develop into quantized vortices in the energy-conserving simulation. 

Next, we consider the dissipative dynamics of the rotating BEC, extending our previous work This paper 

describes more detailed dynamics of the vortex lattice formation by following the time development of the condensate 
density and phase simultaneously. The GP equation with phenomenological dissipation explains the experimental 
results very well. The quasiparticle projection method is also available in the analysis of the dissipative dynamics, 
thus we can study what modes are excited during the dynamical process of vortex lattice formation. In the presence 
of dissipation, the generation of surface ripple is caused by the surface modes with negative frequencies in the rotating 
frame ||]; the onset of this instability is given by the Landau criterion applied to the rotating BEC [J?], ^, |l5[ |l6| . 
Numerical simulation shows that the surface ripples induced by the Landau instability certainly develops to quantized 
vortices. Therefore, it is concluded that vortex nucleation is essentially caused by the Landau instability. The ENS 
experiments jl| should be explained by the two-stage process: vortex nucleation by the Landau instability after the 
thermal component creation by the dynamical instability. 

This paper is organized as follows. Section || introduces our formulation that describes the dynamics of a rotating 
BEC in a harmonic trap. In Sec. OLD, we study the energy-conserving dynamics of a rotating BEC. The dynamics of 



the mode-coupling is analyzed by the quasi-particle projection method. In Sec. IV we study the dissipative dynamics 
of the vortex generation and lattice formation in detail, and make some comments on the origin of the dissipation. 
Our results are compared with the experimental ones. Section |v| is devoted to the conclusion. 

II. THE MODEL 

A. Formulation of the problem 

A BEC trapped in an external potential is described by a "macroscopic wave function" ^(r, t) obeying the GP 
equation. In the frame rotating with the frequency £1 around the z axis the GP equation reads 

lh ^t = (-^ v2 + yt -P + 1/ ™ t -^ + 5l*| 2 -^)*- (!) 

Here g = Air1i 2 a/m represents the strength of interactions characterized by the s wave scattering length a > 0, fi 
the chemical potential, L z = —ih(xd y — yd x ) the angular momentum. The wave function is normalized by the total 
particle number N as J"dr|\l/| 2 — N . An external harmonic trapping potential has the form 

T4 rap (r) = \m{iol{x 2 + y 2 ) + uj 2 z z 2 }, (2) 



and a rotating potential 



Kot(r) = ]^m,uj\{e x x 2 + e y y 2 ) (3) 



with the anisotropy parameters e x ^= e y ; this form describes approximately the rotating potential used in the ENS 
experiments [[jj ^ . Such a rotating potential breaks the rotational symmetry, thus transferring the angular momentum 
into the condensate through the excitation of surface modes or the generation of vortices. 

In order to reduce the system into the two-dimensional x — y space, we separate the degrees of freedom of the wave 
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function as ^(r, t) — ip(x,y,t)<p(z), obtaining the two-dimensional GP equation 



dt 
1 



h 



1 

2m \dx 2 dy 2 



+~mwi{(l + e x )x 2 + (1 + € y )y 2 } 

-M + 9vH(x,y,t)\ 2 - ni 

where 



Hx,y,t), (4) 



and /i includes a constant arising from the integral of 4>(z) . The normalization of the two-dimensional wave function 
ip(x,y) is taken by the particle number N2B per unit length along the z-axis as 

dxdy\^{x,y)\ 2 = N ( [ dz\<t>(z)A = N 2B . (6) 



It is convenient to introduce the scales characterizing the trapping potential; the length, time, wave function are 
scaled as 

t I li> 

x = a h x, t= , i/ , = V^ V 2D — , 

lu± a h 



respectively, with ah = y^h/2mujj_. Then the GP equation is reduced to a dimcnsionless form as 



<9V 



\-(-- 


d 2 \ 


\dx 2 


dy 2 ) 



\{(l + e x )x 2 



-(1 + e y )y 2 } - ii + C|^f - QL 



V>, (7) 



where C = 8irar]N2D and the tilde is omitted for simplicity. 

The two-dimensional approximation may be valid for the condensate in a "pancake-shaped" potential (A = uj z /uj± ^> 
1) or the central part of the condensate in a "cigar-shaped" potential (A < 1). These two types of situations yield 
different forms of the mean field interaction strength C . For A > 1 and huj z larger than the interaction energy, <fi{z) 
is approximated by the one particle ground state wave function in a harmonic potential: 

with atiz = \/h/2muj z = ah/\f\. Then, N213 = N and the parameter C becomes 

C = 8iraN<q = 4%/^AiV— . (9) 

Oh 

On the other hand, for a cigar-shaped condensate with A < 1 one can approximate the system with cylindrical config- 
uration, i.e., translation symmetry along the z direction. By neglecting the spatial derivative term of z component in 
Eq. (|l|) , and the third term in Eq. (||) , the two-dimensional GP equation of Eq. (^) is obtained. Then the parameter 
C is written by 

fWz^dz 8iraN , . 

C = 8 ™^ 2D = ^N { ^ ?dz)2 * _. (10) 



Here R z is assumed to be the Thomas-Fermi radius y / 2fi/muj 2 along the z axis with the chemical potential evaluated 
at the parameter e x _ y — and Q = 0: 

1 r \ 2/5 

15 a x 



H = hw ± [ —NX— ) . (11) 
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The latter approximation is suitable for the ENS experiment |jj which was made under the cigar-shaped potential, 
where C (Eq. (|l0|)) takes the value between 200 and 500. For the large condensate in the MIT experiment 
C r-j 10000, though they did not use the cigar-shaped potential. 

In the two-dimensional analysis, the effect of vortex bending pi] , is not taken into account. Recent experiments 
p6| showed that the time scale of the vortex bending is found to be longer than > 1 sec, which is much longer than 
the time scale of the dynamics of vortex lattice formation (~ 100 msec). We therefore consider our two-dimensional 
analysis effective for the present problem. 

B. Numerical scheme 

The numerical calculations of Eq. (0) are done using an alternating direction implicit (ADI) method [^7j . Defining 
a time step St and space meshes 5 X = 5 y = 5, and denoting the discrete wave function as ip" k — ip(J8 x , kS y , nSt) , ip" k 

develops into V'?^ 1 y i a the intermediate state ip"t 2 as 



and 



+{V jt k + Cm k \ 2 )^ tk } (12) 



At 



T {(^, fe +q^rr)^r 

+ (V jtk + C\^\ 2 )^p}. (13) 

Here we denote At = 6 t /2i, d x ^ k = (^ +1>fc - ^_ 1)fc )/2, d^ k = ^ +1 , fc - 2^ fe +^"_ u and Vj, k = {(l + e x )( 3 5) 2 + 
(1 + e y )(kS) 2 }/4. We used [-128 < j, k < +128] discretized space for the two-dimensional numerical simulation. The 
time step i5 t = 1.0 x 10~ 3 ensures the numerical stability over sufficiently long propagation. 

III. QUADRUPOLE OSCILLATION OF A ROTATING BEC 
A. Time development of the deformation parameter 

We start by discussing the time evolution from the stationary solution in a non-rotating trap. We turn on a rotating 
drive following the experimental procedure of Madison et al. p). The rotation with a frequency VI starts at t = 0, 
and the trap anisotoropy e = {(1 + e x ) — (1 + e y )}/{(l + £x) + (1 + e y)} is increased rapidly from zero to its final 
value 0.025 in 20 msec. The strength of interaction C is set to be 5 00, corres ponding to a = 5.77nm, N = 3 x 10 5 , 
oj z = 11.8 x 27r, A = u)±/lj z = 9.2 ||. The unit for length is — \JU/2moj±_ — 0.728/im and the period of the trap 
9.21 msec. 

Rapid modulation of the trapping anisotropy induces the elliptic oscillation of the condensate. The elliptic oscillation 
is characterized by the deformation parameter ^8| 

a = -n <X l > - <y l > , (14) 

where < A > means J dxdy-ip* Aip . The time evolution of a for several values of £1 is shown in Fig. [I]. For relatively 
small values of O, a makes a simple periodic oscillation with positive values; the initial axisymmetric condensate 
is elongated along the y axis because of the small trap anisotropy e (e x > e y ). As f2 increases from zero, both 
the amplitude and the period increase gradually, leading to the large amplitude oscillation near £1 = 0.7ojj_. For 
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Q = 0.75a->^, however, the periodicity of the oscillation is broken as shown in Fig. [j]. As fl increases further, the 
periodicity is again recovered, but the sign of a changes to negative and its absolute value become smaller and smaller. 
The negative a means that the longer axis of the condensate ellipse is perpendicular to the longer axis of the trapping 
anisotropy. 

This shape oscillation mainly consists of the collective surface mode with angular momentum / = 2, i.e., quadrupole 
mode. In the Thomas- Fermi limit, the dispersion relation for the surface mode reduces to ui = y/lu>± J29|. The 
centrifugal term — Q,L Z lifts the surface mode frequency just as — IVl. For I = 2, hence, it is expected that the 
quadrupole mode is resonantly excited at fl/ujj_ = \/2/2 ~ 0.707. As discussed later, the complete resonance does 
not occur because the quadrupole mode couples with various higher-energy modes through the nonlinear interaction, 
giving rise to the complicated dynamics. We find that for the rotation range 0.72 < Vt/uj± < 0.78 the oscillation 
becomes irregular, reflecting the nonlinear character of the dynamics. The deviation from the pure resonance frequency 
fl/u>± = 0.707 is due to the effect of the trap anisotropy and the nonlinear interaction (see Eq. (|2^)). 

Figure[| shows the profile of the condensate density \ip(x, y, t)\ 2 and the phase 9(x, y, t) = tan -1 (Inx!/>/ReV>) (0 < 6 < 
2ir) when the irregular oscillation occurs. On the surface of the condensate, there appear surface ripples, that is violent 
density fluctuations with a short wave length. In addition, the phase profile shows that many phase singularities, i.e., 
quantized vortices, come into the condensate surface, while the phase keeps the form of the quadrupoler flow inside 
the condensate || psfl . Since these singularities are on the outskirts of the condensate where the amplitude \ip\ is very 
small, they hardly contribute to both the energy and the angular momentum. Such phase singularities, named "ghost 
vortices" in our previous paper |Iif| , are located in the density hollows produced by the surface ripples. 

In the energy-conserving dynamics, such an irregular dynamics is caused by the dynamical instability associated 
with the imaginary frequency of the excitation modes. Sinha and Castin ]l2] | made a linear stability analysis of an 
oscillating condensate with a quadrupolar ansatz, and found the growth of the fluctuation written by polynomials of 
degree n = 3 around Q/loj_ = 0.73. They proposed that the associated dynamical instability triggers vortex nucleation 
in a condensate rotated in an anisotropic harmonic potential, in excellent agreement with the experimental results 
of ENS 0. However, Fig. 2 suggests that the generated surface ripples should be described by polynomials with 
higher-order degrees than n = 3. Thus, we also make a linear stability analysis for n — 4 and 8 by following Ref. 
p2[ , and find that the growth rate of these higher order modes is as much as n = 3. Hence, the dynamical instability 
excites such higher order excitation mode, generating the surface ripples. These ripples slightly increase of the total 
angular momentum because of the presence of ghost vortices, but they never penetrate the inside of the condensate 
to form a lattice. Therefore, it is concluded that the dynamical instability does not lead to "vortex nucleation" as 
observed in the ENS experiments. In Sec. |y|, we will show that the dissipation-assisted instability can make the 
surface ripples develop into the quantized vortices. 



B. Quasiparticle projection method 

The dynamics of the elliptic oscillation is well understood by decomposing the whole dynamics into an assembly 
of fundamental excitation modes. The quasiparticle projection method, developed by Morgan et al. to study the 
nonlinear mixing of the collective excitations ^2|, enables us to decompose a wave function into a condensate and 
noncondensate modes, and monitor the time evolution of their populations. Here we will use this method to study 
the time development of the surface modes excited by the anisotropic rotating trap. 

To construct the mode functions for the projection, we use the solution of the time-independent GP equation with 
the non-rotating axisymmetric trap 



dr 2 



ld_ 

r dr 



r 

T 



i> g {r) = ni> g {r), 



(15) 



where ip g corresponds to the initial non- vortex state in our simulation. The quasiparticle mode functions Ui(r) 
Ui(r)e and Vi(r) = Vi(r)e are obtained by the Bogoliubov-de Gennes equations 



dr 2 



1 d 

r dr 



I 2 



r2 +--p + 2C\4, g \< 



ld_ I 2 

r dr v 



Ui(r) 



+C^ 2 g Vi(r) =uj t u t (r), 



Vi(r) 



+C'i>* 2 Ui(r) = -UiVi(r). 



(16a) 



(16b) 
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The mode functions are subject to the orthogonality and symmetry relations 

d 2 r{ Ul {r)u*{v) - Uj(r)«;(r)} = (17a) 
d 2 r{ui(r)v*-(r) - Vi(r)u*(r)} = 0. (17b) 



Following the method of Ref . ]2^] , we introduce a set of excitations which are orthogonal to ip g . This is achieved by 
projecting out the overlap with ip g from the solutions of the Bogoliubov-de Genne equations, the modified quasiparticle 
wave functions being defined by 

Ui(r) = Ui(r) - aip g {r) (18a) 
V*(r) = <(r) +C *V s (r), (18b) 

where Ci = J d 2 v \_4>gUi\ — — J d 2 r [ipgVi]. The orthogonal relations Eq. ( |l7| ) still hold for these modified wave 
functions. The wave function can be expanded as 

<P(r,t) = {1 + b g (t)}i> g (r) + {Ui(r)bi(t) + v*(r)b*(t)} . (19) 

i>0 

It is easy to show that 1 + b g and bi satisfy the relations 

1 + &,(<) = J d 2 r ^(rty(r,t), (20) 

h(t) = J d 2 r{u*(r)<P(r,t)-v*(r)r(r,t)}. (21) 

The populations of the ground state and the excitation are given by |1 + b g \ 2 and \bi\ 2 . 

In the following discussion, we take only those mode functions that carry angular momentum I but do not possess 
any radial node, i.e., surface mode. Here the index i of Ui and Vi is replaced with I. For I ^ the overlap integral 
ci vanishes because d9e ue = 0, so that ui — ui and vi = vi. Since our system has even-parity for the spatial 
coordinate, no surface modes with odd I are excited. 



C. Fermi- Pasta-Ulam recurrence and chaotic dynamics 

By using Eq. (|Tj), we project the time evolution of the excitation modes \bi\ 2 with I = 2,4, • • -,20 from ip(r,t). 
Figure |^(a) shows the time evolution of I&2I 2 , | £>4 1 2 , |&6| 2 and |1 + b g \ 2 for = 0.7ujj_; other \b{\ 2 are extremely small. 
The rotating potential of Eq. (||) excites most intensively the / = 2 mode which causes the elliptic deformation. The 
small population of the higher modes (I = 4,6, • • •) also appear following the increase of | £» 2 1 2 . Then populations of 
those excited modes return to the initial values almost completely. This simple recurrence repeats periodically in 
time. 

For £1 = 0.75w^ the evolution is somewhat complicated as illustrated in Fig. ||(b). The complete recurrence is 
broken, but some quasi-periodicity still remains. As compared with the simple recurrence dynamics, the depletion 
of the ground state population is remarkable, which means that various modes with higher energy are more excited. 
The superposition of such higher-energy modes produces the surface ripples in the condensate density as shown in 

In general points of view, we face the problem known as Fermi, Pasta and Ulam recurrence phenomenon ||23| . They 
studied statistical behavior in the chain of nonlinear ly coupled oscillators, and found a quasiperiodic behavior of this 
system characterized by returns of the energy to the initial excited mode. Later, it was shown that there exists a 
threshold for the onset of "stochastization" which is brought by high-energy excitations. For BEC in an anisotropic 
potential, the nonlinearities inside the condensate may give rise to the stochasticity in its time evolution J|(| , which 
has been found in the numerical simulation of the GP equation |3l], ^] . 



D. Two- mode analysis 



The simple analysis of the equation of motion for the quasiparticle population helps us understand the behavior of 
Fig. H (a) and (b). Substituting Eq. (|l|) into Eq. (Q), we get 
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.db a 



dt 



^ = / d 2 rr g V m t{% + A) + C / d 2 r 



(fi B + W„| 4 + 2| A^| 2 + A^f + A|A| 2 ^ 



Mi 

! 9F 



(w, - ifi)6i(*) + / d 2 r<y rot (^ g + A) + C / d 2 ™; 



2|A| 2 ^ + A 2 C + A|A| 2 + |V fl | V 9 (& s + K) 



+ / d 2 n;f Kot(^ 5 + A*)+C d 2 rv\ 



2|A|>* + A*> ff + A*|A|^ + M 2 r Q (b g + b 



(22) 



(23) 



where A = bgtp g + J2i{ u i( r )bi(t) + v i ( r )fy* (t)} and V rot = (e x x 2 + e y y 2 )/4:. The presence of V T0 t makes the integrals 
/ d 2 rip* V ro tUi or J d 2 rui,V ro tVi ■ ■ ■ finite, and these terms are reduced to the form To^z,/' + T\8ij,i-i + T 2 6u>+ 2 after 
the integral of ^-component (if e x = e y , the off-diagonal terms vanishes). Hence, the ground state is coupled directly 
to only the I = 2 mode in Eq. ( ^2|) through V ro t- The remaining terms with C couple the various modes with the 
same parity to each other, resulting in a complex time evolution. 

As seen in Fig. |](a), the dominant contribution to the dynamics comes from the change of the population |1 + b g \ 2 
and |&2| 2 - To understand the basic properties we take only terms with I = and 2 (two- mode approximation). In 
addition, we neglect the contribution of i>2(f*), because this term represents the excitation traveling oppositely to the 
rotation. Then Eqs. (p2]) and (23) reduce to 



■db g 
.db 2 



Tbo(l 



P 00 {2Re(bg) + \bg\ 2 }(l + bg) 



+Po2(l + b g )\b 2 \ 2 + T 02 b 2 , 
{lu 2 - 2n + T 22 )b 2 + P 22 \b 2 \ 2 b 2 
+Po 2 {2Re(6 5 ) + \b g \ 2 }b 2 + T 02 (l + b 



(24a) 
(24b) 



where 



d 2 rip*V xo tU2 = I d 2 ru* 2 V mt il)g, 



Too = I d ri(i*V Iot ipg, T 22 = \ d ru* 2 V rot u 2 
Poo 

Pq2 



C J d 2 r\^ g \\ P 22 = C I d 2 r\u 2 \\ 
2C J d 2 r\4, g \ 2 \u 2 \ 2 . 



It is convenient to represent the complex values b g and b 2 in terms of the amplitude and the phase as = |1- 



and b 2 = \b 2 \e 2 . Since the total population |1 



population difference p 



1 



b g \ 2 



bg? 



-bgW 



\b 2 \ is a constant of motion, Eqs. (E4j) have two variables: 



1 62 1 an d relative phase 6 



Then, Eqs. (M reduce to 



dp 
~db 



d6 



2p 



p 2 sin ( 



— = Up- T 02 r 

dt \/\-p 2 



cos( 



Au 



(25) 
(26) 



with the conserved Hamiltonian 



H(p,e) = -u P 2 



2T 02 \J\-p 2 cos 9 - Aojp, 



(27) 



where U = (Poo + P22 — 2Po 2 )/2 and Aui — lo 2 — 2Vl + T 22 — Too + U. These formula are the same with those used 
in the Josephson dynamics of the condensate in a double well potential, and the exact solution of p(t) is expressed 
by the Weierstrassian elliptic function (3^]. For C = 500, we numerically determine as lo 2 — 1.438a; j_, U = —0.05, 
Slu = cu 2 — 20 + 3e — 0.05 and T 2 = 1.14e from the calculated ip g and u 2 . The solution of Eqs. fl25| ) and ( p6| ) with 
these parameter values are shown in Figs. ||(c). The periodicity and the amplitude of |1 + b g \ 2 reproduce well the 
results of the GP equation. 

Figure |] shows the dependence of the maximum of the \b 2 \ 2 oscillation on Ct/u>j_- The maximum grows near 
fl = 0.7uj±. Note that as e increases the peak is shifted to the larger value of fi - from the pure I = 2 resonance 
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frequency £1 = = 0.719ui±. The dynamics of the GP equation does not cause pronounced resonance of I62I 2 

because the transition into the higher energy modes occurs. For e = 0.025 we found that the chaotic oscillation occurs 
in 0.72 < il < 0.78, which corresponds to the range where the maximum of I&2I 2 is more than 0.4 in Fig. |[ Following 
this criterion, we obtain the £1 — e parameter region where the large amplitude oscillation is expected, as shown in 
the inset of Fig. ||. In this region the mode coupling via the nonlinear interaction becomes important, leading to a 
chaotic dynamics. We confirm by the numerical simulation that the simple recurrence is indeed broken in this region. 
Such a nonlinear mode c oupli ng seems to be important in the experiment on vortex nucleation W , and the detailed 



will be discussed in Sec. IV E 



IV. DYNAMICS OF VORTEX LATTICE FORMATION 

In this section we focus on the dissipative dynamics of a trapped condensate following the sudden turn-on of rotation. 
The dissipation is necessary to simulate the dynamics of vortex lattice formation by the GP equation, because a vortex 
lattice corresponds to a local minimum of the total energy in the configuration space [ p4| [35| . Preliminary results 
were reported in Ref. [ fl4"| . This section deals with the detailed dynamics by following the time development of the 
condensate density and the phase simultaneously, and presents a number of previously unpublished results. Although 
the dissipation is treated phcnomcnologically in the GP equation, the simulation explains the experimental results 
very well. 

The results show that the generation of surface ripples is also achieved by the instability of negative eigenvalue 
modes, i.e., Landau instability, and the following time development certainly leads to vortex penetration into the 
condensate. The experimental results on vortex nucleation by the ENS group jj], ^ can be understood by taking into 



account both the dynamical instability and the Landau instability. This is described in Sec. IV E 



A. Phenomenological dissipative equation 

Before discussing the detailed dynamics, we make some comments on the dissipation. As in the previous study E3, 
the dissipation is treated phenomenologically in the GP equation. The time derivative term of Eq. (0) is modified as 



L (*__ 


d 2 \ 


\dx 2 


[ dy 2 ) 



+ {l + e y )y 2 }- n + C\ij\ 2 -m z 



V>, (28) 



where the dimensionless parameter 7 describes the dissipation. This form of the dissipative equation follows the study 
of Choi et al. and that of other superfluid systems |37j]. Choi et al. determined the value of 7 to be 0.03 by 
fitting their theoretical results with the MIT experiments on collective damped oscillations |3^| . Thus, we also use 
7 = 0.03 throughout this work. Since this dissipative term is much smaller than other terms in the GP equation, a 
small variation of 7 does not change the dynamics qualitatively but only modifies the relaxation time scale. 

The phenomenological dissipative equation ( p8| ) may be relevant to the recent numerical work by Jackson and 
Zaremba [ |39| on the coupled dynamics of a condensate and a noncondensate. Their simulation is based on the the 
generalized GP equation at a finite temperature 

( h 2 v 2 \ 
ih—=l-—— + V tIap +gn + 2gh-iTW. (29) 

This equation was derived by Zaremba, Nikuni and Griffin j40|] , where n(r,t) = \^(r, t)\ 2 is the condensate density 
and n(r, t) the noncondensate density. The dynamics of the noncondensate was described by the Boltzmann kinetic 
equation for the distribution function of the noncondensate atoms. The noncondensate atoms are assumed to obey 
the single-particle Hartree-Fock spectrum in their formulation. The numerical simulation by Jackson and Zaremba 
|f39f well explains the experimental results by the JILA group jll| . 



Equation (|28j) can be derived from Eq. (29) with some additional approximations. We treat the noncondensate as be- 
ing in static thermal equilibrium, and neglect the mean field of the noncondensate under the assumption n n. Then, 
the dissipation of the condensate motion is associated with the term T(r,t) = (h/2n)Ti2 — (%/2n) f{dp/(2irTi) 3 }Ci2 
with the collision integral Cyi between the condensate and noncondensate atoms. Under the local equilibrium distri- 
bution of the thermal atoms, Y\i is proportional to the difference of the local chemical potential between condensate 
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and noncondenate as ]?i2 oc fi nc (r, t) — f-i(r,t) pQ|. Approximating (i(r, t)ty ~ —ih(d^/dt) p3[ , we obtain 

(i - 7 )ft-^ = (^-^V 2 + y trap + 5 |*| 2 - i7 MncJ *• (30) 

Compared with Eq. (|2^), the chemical potential of this equation is replaced by that of noncondensate fi nc , yielding 
the imaginary term i 7 /Lt nc . Note that the relation fi — /i nc is satisfied for the equilibrium condensate. If the space 
and time dependence of ^ nc is neglected, this equation becomes Eq. ( p8| ) by the transformation ^ — > vfe" 4 ^'/^ and 

Mnc -> 

The time development of Eq. J2g|) does not conserve the norm of the wave function as well as the energy. In our 
simulation, the chemical potential fi is adjusted at each time step in order to conserve the norm and to decrease 
the total energy monotonically. Recently, Eq. ([30j) was also derived by Gardiner et al. by another approach [Elj, 
extended to the simulation of vortex lattice formation from the rotating thermal cloud iffBf . They made the numerical 
simulation with the fixed chemical potential, finding that the norm of \& decreased or increased as the system evolved. 
In the actual experiment, such a loss of the number of condensate atoms is caused by the rotation- induced heating [0, 
so that the generated thermal component might affect the dissipative dynamics of the vortex generation. However, 
under the phenomenological model, the change of the norm and the time scale of the dynamics is only a few % by the 
effect, whether the chemical potential is fixed or not. The detailed quantitative study of those processes is beyond the 
scope of this paper; our treatment for the chemical potential is adequate to describe the actual dissipative dynamics. 

The assumption of the static noncondensate may be applicable to the experimental condition in Ref. & ^J, 
According to the estimation by Guery-Odelin [E6[ , the spin up time for the whole noncondensed atoms to catch up 
with the rotating trap is about 15 seconds in collisionless regime. Since the typical time for the vortex formation of 
the condensate is a few hundred milliseconds, the condensate motion is separated from the noncondensate one under 
the rotating perturbation. However, the dynamic coupling of the mean-field between condensate and noncondenate 
causes the condensate motion to be damped (known as Landau damping), which is not included in Eq. (p8|). Such a 
damping is several times larger than that by the Ci2-collision at very low temperature |3^, ^7|. Thus an additional 
damping may be provided by the full coupled dynamics of a condensate and a noncondensate. 

The value of 7 is estimated by following the formulation for a uniform Bose gas p2| . The parameter 7, in which 
length and energy are scaled by ah and huj±, has the form 



'■^(tJVW (31) 

Here A is a factor of order unity for T > 0.5Tc and approaches zero as T — > 0. Using the typical experimental 
parameters, for example, T = O.bTc, Tc = 500nK, a = 5.5nm, n ~ 10 14 /cm 3 and oj± = 100 x 27rHz, we obtain 
7 ~ A x 10~ 2 which is consistent with 0.03 used in this paper. 



B. Dynamics from a non- vortex state to a vortex state 

Using Eq. (^8|), we discuss the dynamics of vortex lattice formation in more details. As in the previous section, a 
sudden switch-on of rotation is made for the condensate with C — 500. Figure || shows the time development of the 
condensate density \ip(x, y, t)\ 2 for Q/u>± = 0.7 |fj8| . Initially, the condensate makes a quadrupole oscillation, but the 
oscillation is damped due to the dissipation. After a few milliseconds, the boundary surface of the condensate becomes 
unstable, generating the surface ripples which propagate along the surface as shown in Fig. |^(c). The excitations are 
likely to occur on the surface whose curvature is low, i.e. parallel to the longer axis of the ellipse. Then the waves on 
the surface develop into the vortex cores in a very short time [Fig. ||(d) and ||(e)]. As is well known in the study of 
rotating superfluid helium p4 
between vortices tends to pus 



the rotating drive pulls vortices into the rotation axis, while repulsive interaction 
them apart; this competition yields a vortex lattice whose vortex density depends on 
the rotation frequency. In the presence of dissipation, six vortices enter the condensate, eventually forming a vortex 
lattice. As the vortex lattice is being formed, the axisymmetry of the condensate shape is recovered. 

The rotating potential V ro t has even-parity with respect to the coordinate. Accordingly, the dynamics is symmetric 
in the x-y plane, and the number of the generated vortices is always even. The even-parity forbids the dynamics in 
which the system develops an asymmetric steady state. To remove this restriction, we introduced an infinitesimal 
artificial perturbation with odd-parity in Vtrap- This perturbation allows the system to develop into an asymmetric 
steady state as shown in Fig. ||(h), the energy of which is lower than that in Fig. ||(g). It takes a few hundred msec 
for the transition from Fig |5|(g) to ||(h) . 

The corresponding time development of the phase of ip{ x i Vi t) is shown in Fig. |^. As seen in the energy-conserving 
dynamics, as soon as the rotation starts, the phase field inside the condensate takes the form of quadrupolar flow 
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9(x 1 y) = axy + const, and just outside the Thomas-Fermi boundary there appear ghost vortices; for example, Fig. 
||(b) shows about 20 vortices. Ghost vortices move toward the rotation axis, but their invasion into the condensate is 
prevented at the Thomas-Fermi boundary. However, as the surface ripples are generated, the ghost vortices start to 
penetrate the condensate. There takes place the selection of the defects to penetrate, because their further invasion 
costs energy and angular momentum. For example, Fig. ^ shows that six vortices enter the condensate and form a 
lattice, while other excessive vortices are repelled and escape to the outside. 

As seen from Fig. |j| vortex invasion is likely to occur on the surface parallel to the longer axis of the ellipse. This 
is simply understood by the velocity field of the elliptic condensate v = Waxy — (ay, ax) as seen in Fig. ||(b). Near 
the condensate surface parallel to the longer axis of an ellipse, this velocity field has the same direction as the velocity 
field made by the ghost vortices. There the additive velocity field v works as the attractive force which pulls the ghost 
vortices into the condensate. While the velocities of the condensate and the ghost vortices have opposite direction 
near the surface parallel to a shorter axis, where the condensate dislikes the invasion of the ghost vortices. Therefore, 
the vortices enter the condensate from the surface parallel to the longer axis. 

Figure [7] shows the time evolution of the deformation parameter and that of the angular momentum per atom 
tzj^i — J dxdyip*(L z /fi)ip in our dynamics of Fig. ^| and [| They very well reproduce the experimental results of 
Ref. @. Before 300 msec, both a and £ Z /Ti make damped oscillations. When vortices enter the condensate, a falls 
abruptly to a value below 0.05 and £ z /h increases to 4 reflecting the number of the generated vortices. 

The final equilibrium value of £ z depends on the number of the vortices which form a lattice. Figure || shows the 
dependence of the number of vortices on the frequency fl/to± and the angular momentum per atom £ Z /Ti for C =250, 
500 and 1500 at 800 ms after the rotation starts. The increase in C stabilizes the lattices of more vortices for the 
same frequency, and reduces the critical frequency at which the first vortex appears. Note that the value of £ z /h is 
about a half of the number of vortices. This is understood by the simple model in which the condensate with a vortex 
lattice makes a rigid-body rotation. Then the mean angular momentum per atom at r — \/ x 1 + y 1 is £ z /h — mflr 2 . 
The average of the angular momentum per atom averaged over the whole condensate is given by 

<//ft>-/JM* (32) 

Assuming the spatially homogeneous density, we obtain < £ z /fi >— mflR 2 /2h with the typical radius R of the 
condensate. In the limit of a rigid-body rotation, the number of vortices i\^, attlco at the rotation frequency f2 is given 
by Feynman's rule 

iV-J-ttice = ^ = mfi# 

k n 

where n v represents the number of vortices per unit area and k the quanta of circulation. Hence, one obtains 
< £ z /h >= Ay attlcc /2. The numerical result better agrees with this estimation for larger and larger C because the 
condensate with a dense vortex lattice mimics a rigid-body rotation p9| . The small disagreement may be attributed 
to a small deviation from Feynman's rule and the inhomogeneous density. 



C. Dynamics starting from an initial state with one vortex: metastable state 

Next, we discuss the time evolution starting from a one-vortex state. As seen from Fig. [s| and ^, as soon as the 
rotation is turned on to the irrotational condensate, the condensate makes the quadrupole deformation and its phase 
takes the form 9(x, y) — axy + const. This behavior will be changed if the dynamics starts from the initial state 
with one vortex which has already the circulating phase field. References |ll| and discuss the stability of the 
condensate with a vortex for the quadrupole mode, and connect it with the additional vortex formation. To investigate 
this problem, we prepare the initial state with one vortex for C = 500 and = 0.4wj_ larger than the thermodynamical 
critical frequency stabilizing one vortex state || S, 11, and start to rotate the system with fl = 0.7uj± as before. 
The time evolution of the phase is shown in Fig. |)|. The numerical simulation reveals the nontrivial structure of the 
phase field; at the center of the condensate the phase maintains the circulation carried by an original vortex, while 
in the outer region the phase makes the quadrupolar flow [Fig. |9|(b)]. Therefore, the condensate with one vortex also 
makes quadrupole deformation. The corresponding time evolution of the deformation parameter a is shown in Fig. 
|?] by the dotted line. The small amplitude of a compared to the previous result is due to the shift of the resonance 
frequency of the quadrupole mode because of the presence of a vortex Q . 

After that, the dynamics follows the same process as before. The final steady state is the lattice with seven vortices. 
This state is energetically higher than that of Fig. |5|(h) with six vortices. Therefore, by starting from different initial 
states, one can obtain various metastable states with different configuration of a vortex lattice, as studied in rotating 
superfluid helium pij. These metastable states are observable experimentally. 
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D. Surface ripple excitation via the Landau instability 



As stated above, vortices are initially generated outside the Thomas-Fermi boundary of the condensate, where 
the energy cost to create defects is small because of the extremely low density. However, their penetration into 
the condensate was accomplished with the help of the surface ripples, induced by the instabilities in the non-vortex 
state. The quasiparticle projection method is useful to reveal the instability of the surface mode in this dissipative 
dynamics. Note that the time derivative term of Eq. ( p2| ) and (53) is modified as (i — j)d/dt. Then, the non- vortex 
state becomes unstable when at least one excitation frequency H>i = u>i — Id + 0(e) becomes negative, causing the 
exponential growth like bi(t) ~ e~ 7 " ! *. Isoshima and Machida examined that the instability associated with the 
negative excitation frequency gives rise to the vortex formation |8| , and calculated the critical frequency O c at which 
the first vortex appears within the Bogoliubov theory. Garcia-Ripoll and Perez-Garcia calculated £l c for more realistic 
conditions (ll[] . That critical frequency can be expressed by the Landau criterion applied to the rotating BEC [Q : 

n c = minf^Y (34) 



v I , 

The angular momentum l c which yields Q c takes a value more than 4 with the parameter used in experiments [f^, |l5| ; 
for C — 500 used in this paper, l c =8 and fi c = 0.5lu±. 

Our simulation confirms that this instability actually leads to the vortex generation, where the hollows of the surface 
ripples always evolve into the vortex cores. The vortex core near the condensate surface has the size of the coherence 
length £ determined by the local density. The number of the vortices generated at the condensate surface may be 
approximately given by N* u ~ 2ttR/^. The numerical solution shows 2ttR/^ is nearly equal to l c . This is understood 
by the fact that the surface mode with l c has the wavelength of the order of £ as discussed in Ref. . 

Note that N^ mi differs generally from the number of vortices jV^ attlce , depending on fi, in an eventual vortex 
lattice. This fact classifies the dynamics of a vortex invasion into two regimes. When > 7V^ attlco , the vortices 

which invade the condensate are chosen from iV™ rf vortices generated at the surface and form a lattice following the 
dissipative vortex dynamics, the extra N^ ui:l — ]\[^ ttlcc vortices being expelled out j35|. This dynamics is shown in 
Fig. |l^(a) in terms of the quasiparticle populations for N c ~ l c = 8 and 7\Tj, attlcc = 2; for £1 = 0.57lu± the excitation 
frequencies lui with I — 4 ~ 14 are negative. At the moment the vortices are about to enter the condensate (t ~ 1.0 
sec), the modes with I =4, 6 and 8 (with \bi\ 2 > 0.02) are excited. These modes rapidly decay as two vortices penetrate 
into the condensate. After the vortex invasion, the final growth of the projected population of the / = 2 mode reflects 
the phase structure in the surface region of the lattice with two vortices. 

On the other hand, if N* ulf < j\r|, attlce , the number of the vortices first generated at the condensate surface is not 
sufficient to form a final lattice, so that successive invasion of vortices is needed. This situation is shown in Fig. |lO|(b) 
with jV,Jf ttlcc = 14. The frequency a>2 is already negative at f2 = 0.86wj_, which causes the exponential growth of ^2 1 2 ■ 
When the condensate deforms elliptically at t ~ 0.2sec, many kinds of surface modes are excited violently. Then 
the quasiparticle populations with high angular momentum (I = 10,12 and 14 are plotted in Fig. |l^(b)) oscillates 
during a short time after this burst. Finally |&i4| 2 grows when 14 vortices enter the condensate completely. Note that 
the I = 14 mode does not grow shortly after the burst, but it grows gradually through the excitation of the lower 
angular momentum modes. The observation of the time development of a density field could find that the radius of 
the condensate is increased unsteadily with the successive invasion of vortices. 



E. Relation between dynamical instability and Landau instability in the experiments on vortex nucleation 

For an irrotational BEC subject to rotation, we have clarified by the numerical simulation of the GP equation 
that there exist two instabilities that is relevant for vortex generation. The dynamical instability appears when 
the quadrupole surface mode is resonantly excited. The Landau instability associated with the negative excitation 
frequency is effective only in the presence of the dissipation. In this section, we discuss which instability is important 
for actual vortex lattice formation by referring to the experimental results. 

In the experiments |l], |[ |), the authors observed that the vortices are nucleated most near = 0.7cjj_. Thus, 
vortices are generated when the condensate becomes resonant with the rotating perturbation that excites a quadrupole 
mode. These results strongly support the dynamical instability scenario by Sinha and Castin However, they 

also observed the vortices at lower off-resonant frequency, which cannot understand by the dynamical instability. On 
the other hand, the critical frequency at which the first vortices appear is extensively discussed by several authors, 
based on the Landau criterion. M ||, O, O, |l6|, [l7], [l8|]. The relation between the dynamical instability and Landau 
instability is still controversial, thus we will discuss these relation in the range of on- and off-resonant frequency 
separately. 
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In the on-resonant range 0.72 < fl < 0.78, the surface ripple is excited, but the dissipation is necessary to generate 
vortices in a condensate. Here the dissipation originates in the thermal component, which should be almost negligible 
in the experiment of the atomic gas at very low temperatures. However, in the dynamical process of a condensate, 
there is a possibility that the thermal component will be produced under a strong perturbation. The experimental 
result of ENS m may be explained as follows. Consider a situation in which there is almost no dissipation at very low 
temperatures, namely, Landau instability does not work there. In a quadrupole resonance, however, the dynamical 
instability causes stochasticity in condensate oscillations, leading to the creation of the thermal component |3l], |32| . 
Indeed, Hodby et al. in Oxford has reported that temperature of the system increases from 0.5T C to 0.8T C during the 
vortex formation procedure J3|. The time spent in this process is determined by the growth time of the dynamical 
instability, which is expected to be about 100 msec |Q. Then, the created thermal component makes the dissipation 
effective, vortices penetrating into the condensate via the dissipation-assisted Landau instability. To make clear this 
hypothesis we need the analysis beyond the mean field and leave this issue for future study. 

In the off-resonant range, a condensate makes only a stable quadrupole oscillation without dissipation. However, 
our results show that vortices may be generated beyond Q c whenever finite dissipation works. Therefore, if we make 
the experiment in which the temperature is so high that the dissipation works effectively, it is possible to observe 
the critical frequency given by Landau criterion. For the weakly anisotropic rotating potential used in the ENS 
experiments and in this paper, excitation of high-Z modes can be made only through the I = 2 excitation, so that it 
takes very long time for the vortex formation near the critical frequency. Figure O represents the relaxation time 
for vortex lattice formation in our simulation with C = 500 and 7 = 0.03 (f2 c = 0.5). Each relaxation time is taken 
at the moment when the angular momentum becomes a half of the final equilibrium value during the rapid increase 
(see Fig. ^). We also show the growth time via Laudau instability r ; = —{7(0;; — Ityuj^}^ 1 for I = 2, 4 and 6. For 
51 < 0.7a; j_, where 0J2 — 2f2 > 0, the relaxation time is longer than 77 with / = 4, 6, 8,..., because high-/ mode is only 
excited through the I = 2 mode which has no instability. Near the critical frequency fi c one finds too long relaxation 
time about 1 sec. As fi increases higher than 0.8wj_, where 0^2 — 217 < 0, the relaxation time matches t^. The ENS 
group [j| did not continue the observation beyond 1 sec, finding no vortices in the off-resonant range near fl c . The 
longer rotation and the presence of the dissipation could nucleate vortices in the rotation range beyond Q c . 

V. CONCLUSION 

We investigate the detailed dynamics of a rotating BEC in a trapped condensate following the sudden turn-on of 
rotation. The numerical analysis of the two-dimensional GP equation shows a series of nonlinear dynamics that has 
not been clarified so much. In the energy-conserving dynamics, we study the quadrupole oscillation induced by the 
anisotropic rotating potential, and the time development of excitation modes by the quasiparticle projection method. 
In the resonance range of the quadrupole mode, the large amplitude oscillation causes the nonlinear mode coupling 
towards the higher energy modes, reflecting the dynamical instability. This nonlinear process leads to generate surface 
ripples, but not to nucleate vortices. In the dissipative dynamics, the vortex lattice formation is revealed in more detail. 
The vortex penetration into the condensate is achieved by the surface modes excitation with negative frequencies in 
the rotating frame, so that the critical frequency for vortex generation is determined by the Landau instability. Two 
possible instabilities for vortex generation is discussed by comparison with the experiments. Although the dynamical 
instability may promote vortex formation, such an instability alone cannot explain the experimental results. The 
dynamical instability plays the role of increasing a thermal component, and vortex nucleation and penetration is 
caused by the Landau instability. 

The similar situation with the competition between the dynamical instability and the Landau instability has been 
found in the center-of-mass oscillation of a BEC in a one-dimensional optical lattice |nj . They observed the critical 
superfluid speed above which the dissipative dynamics starts. Its main features, including the parameter range of 
instability, can be explained by the simulation of the dissipationless GP equation |)2j , which means that the observed 
phenomena may be caused by the dynamical instability. However, the following development of a condensate is 
not explained only by the dynamical instability; the dissipation-assisted instability such as the Landau instability is 
necessary to describe the dynamics These problems on the dissipative dynamics of a BEC, as well as the vortex 
lattice formation, offers the testing ground for the analysis beyond the framework of the GP equation. 
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FIG. 1: Time evolution of the deformation parameter a for fl/ivx =0.65, 0.70, 0.75, 0.80, 0.85. 



FIG. 2: The density and phase profile of the condensate with Q/ui± = 0.75 at 105 msec. In (b), the value of the phase changes 
continuously from (black) to 2n (white). There appear some lines where the phase changes discontinuously from black to 
white. These lines correspond to the branch cuts between the phases and 2ir, and their apexes around which the value of the 
phase rotates continuously from to 2n represent phase defects, i.e., quantized vortices. 
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FIG. 3: Time evolution of the surface mode populations \bi\ 2 for (a) Q/lj±_ = 0.7 and (b) Q/wx = 0.75. Figures (c) represents 
the time evolution of |1 + b g \ 2 (dashed curve) and |&2| 2 (solid curve) in the two-mode approximation 

FIG. 4: The dependence of the maximum of | £>2 1 2 on Sl/wi in the two-mode approximation for e =0.01, 0.025, 0.05 and 0.08. 
The inset shows the parameter region in which the large amplitude oscillation with max(|&2| 2 ) > 0.4 occurs. 



FIG. 6: Phase profiles of ip corresponding to the density profiles of Fig. |5j. The value of the phase changes continuously from 
(black) to 2-7T (white). The discontinuous lines between black and white correspond to the branch cut of the complex plane, 
and their edges represent quantized vortices. The unit of length is the same as that of Fig. H. 



FIG. 7: Time evolution of the distortion parameter a (solid curve) and the angular momentum per atom £ z /h (dashed curve) 
in the dynamics of Fig §and|§ The dotted curve shows a for the dynamics starting from the initial state with one vortex. 



FIG. 8: The number of vortices (solid curve) and angular momentum per atom (dotted curve) versus Q/u±, for C =250 (with 
empty circles), 500 (with filled circles) and 1500 (with crosses). 



FIG. 9: Phase profile of the simulation starting from one vortex state. 



FIG. 10: Time evolution of the surface-mode amplitude \bi\ in the dissipative dynamics for (a) Q = 0.57u>± and (b) = 0.86u!± 



under the same condition in Sec. 



IV 



In (b), we plot only modes with 1 = 2, 10, 12, 14, which are important in the discussion. 
The inset shows the density profile of the final steady state. The other inset in (b) shows |feio| 2 , I&12I 2 and |fei4| 2 near t = 0.3 



FIG. 11: The circles represent the relaxation time of vortex lattice formation by the numerical simulation for C = 500 and 
7 = 0.03. Thin solid curves show the 77 = — {7(0^ — IQ.)ujj_}~ 1 for I = 2, 4 and 6 obtained from Eq. ( |l^ ) 



FIG. 5: Time development of the condensate density \tp\ 2 after the trapping potential suddenly begins to rotate at t = with 
n = 0.7wl. 



